{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению. Сессия № 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <center> Индивидуальный проект:\n",
    "## <center> <font color='Chocolate'>Прогнозирование продолжительности поездки нью-йоркского такси</font>\n",
    "Автор материала: ***Бадамшин Булат Ильдарович***\n",
    "\n",
    "e-mail: aem-team@ya.ru\n",
    "\n",
    "slack: Bulat Badamshin\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### План.\n",
    "1. Описание набора данных и решаемая задача.\n",
    "2. Первичный и визуальный анализ данных.\n",
    "3. Построение новых признаков.\n",
    "4. Визуальный анализ всех признаков.\n",
    "5. Корреляционный анализ признаков.\n",
    "6. Выбор метрики.\n",
    "7. Выбор модели.\n",
    "8. Настройка гиперпараметров модели с помошью библиотеки hyperopt.\n",
    "9. Построение валидационных, тестовых и тренировочных кривых.\n",
    "10. Стекинг.\n",
    "11. Выводы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Описание набора данных и решаемая задача"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Задача\n",
    "**Прогнозирование продолжительности поездки такси в Нью-Йорке**\n",
    "- Задача является актуальной для оператора такси. Прогнозирование времени поездки позволит оператору  оптимизировать свою работу и улучшить финансовые показатели.\n",
    "- С точки зрения учебных целей эта работа также представляет большой интерес. В задаче предстоит работать с различными типами данных, в том числе с геораспределенными. Эти данные позволяют генерировать большое клоичество новых признаков. Что выделяет одну из особенностей задачи - возможность построения и исследования большого количества новых признаков.\n",
    "- Кроме того, организаторы допускают использование дополнительных данных. Так, например, предлагается использовать информацию о погоде в Центральном праке за 2016 год. \n",
    "- Объем данных достаточно большой (~1.5 млн записей тренировочной выборки). Что с одной стороны вносит дополнительные требования к вычислительным ресурсам, с другой стороны большой объем данных позволит получить устойчивое решение."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Описание данных"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Данные взяты из kaggle-соревнования [New York City Taxi Trip Duration](https://www.kaggle.com/c/nyc-taxi-trip-duration). \n",
    "\n",
    "[Ссылка на данные](https://yadi.sk/d/TSQV2HeU3Pd56P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Датасет представлен в виде csv-файлов для train и test-выборки. Данные соревнования основаны на источнике о поездках [NYC Yellow Cab](https://cloud.google.com/bigquery/public-data/nyc-tlc-trips) в 2016 году.\n",
    "- Train-выборка включает 1458644 записей.\n",
    "- Test-выборка включает 625134 записей.\n",
    "- Целевая переменная **Trip duration** - продолжительность поездки в секундах.\n",
    "- Признаки:\n",
    "    - **id** - уникальный идентификатор поездки;\n",
    "    - **vendor_id** - идентификатор провайдера;\n",
    "    - **pickup_datetime** - дата и время включения счетчика;\n",
    "    - **dropoff_datetime** - дата и время отключения счетчика;\n",
    "    - **passenger_count** - число пассажиров в машине;\n",
    "    - **pickup_longitude** - долгота точки, где был включен счетчик;\n",
    "    - **pickup_latitude** - широта точки, где был включен счетчик;\n",
    "    - **dropoff_longitude** - долгота точки, где был выключен счетчик;\n",
    "    - **dropoff_latitude** - широта точки, где был выключен счетчик;\n",
    "    - **store_and_fwd_flag** - этот флаг показывает была ли задержка в передачи данных о поездке вендору: Y - задержка передачи данных ведору была, N - задержки не было.\n",
    "\n",
    "В задании подразумевается использование все доступных дополнительных признаков."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 2. Первичный анализ данных "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import re\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
    "from sklearn.model_selection import train_test_split \n",
    "from sklearn.model_selection import cross_val_score, KFold\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.linear_model import SGDRegressor, LassoCV\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn import metrics\n",
    "from sklearn import grid_search\n",
    "from sklearn.cluster import KMeans, MiniBatchKMeans\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "from sklearn.metrics import mean_squared_error, mean_absolute_error\n",
    "import sklearn.metrics\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "import lightgbm as lgb\n",
    "#import xgboost as xgb\n",
    "\n",
    "import pickle\n",
    "import tqdm\n",
    "\n",
    "from collections import Counter\n",
    "\n",
    "import statsmodels.api as sm\n",
    "from scipy.stats import shapiro\n",
    "\n",
    "RANDOM_STATE=0\n",
    "\n",
    "%pylab inline "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Загрузим обучающую и тестовую выборки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train = pd.read_csv('../../data/train_NYC_Taxi.csv')\n",
    "test = pd.read_csv('../../data/test_NYC_Taxi.csv')\n",
    "print (train.shape, test.shape)\n",
    "train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- train-выборка включает в себя 11 признаков, из которых один целевой. \n",
    "- test-выборка включает в себя 9 признаков. Признака dropoff_datetime нет в тестовой выборке."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.info(), test.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train['store_and_fwd_flag'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Обучающая и тестовая выборки не имеют пропусков.\n",
    "- Признак id можно считать категориальным, но скорее всего он не несет полезной информации.\n",
    "- В обучающей выборке 2 признака типа timestamp, в тестовой один.\n",
    "- В обеих выборках присутствует бинарные признаки: store_and_fwd_flag и vendor_id.\n",
    "- Остальные признаки вещественные."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Проанализируем целевую переменную. \n",
    "#### Выбросы по целевой переменной trip_duration\n",
    "- Для анализа выбросов построим вариационный ряд по целевой переменной"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "plt.scatter(range(train.shape[0]), np.sort(train.trip_duration))\n",
    "plt.xlabel('index')\n",
    "plt.ylabel('trip duration')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax =plt.subplots(figsize=(8,6))\n",
    "ax = train['trip_duration'].hist(bins=100)\n",
    "ax.set_title('Распределение целевого признака trip_duration');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- На графике вариационного ряда видно, что присутствуют выбросы - аномальные продолжительности поездки. \n",
    "- Распределение сильно скошено и выглядит как константное из-за больших значений аномалий."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Выкинем из обучающей выборки записи, целевая переменная в которых не попадает в 99% распределения."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = train.trip_duration.quantile(0.99)\n",
    "train = train[train['trip_duration'] < q]\n",
    "\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.scatter(range(train.shape[0]), np.sort(train.trip_duration.values))\n",
    "plt.xlabel('index')\n",
    "plt.ylabel('trip duration')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax =plt.subplots(figsize=(8,6))\n",
    "ax = train['trip_duration'].hist(bins=100)\n",
    "ax.set_title('Распределение целевого признака trip_duration');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- После фильтрации выбросов целевой переменной вариационный ряд представляет собой монотоную линию без разрывов (выбросы устранены).\n",
    "- Распределение стало менее скошенным. Теперь с ним можно работать."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Проверка на нормальность распределения целевой переменной с помощью критерия Шапиро-Уилка и QQ-плота"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tg = train['trip_duration']\n",
    "sm.qqplot(tg)\n",
    "pv = shapiro(tg)[1]\n",
    "plt.title('p-value: ' + str(pv));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Гипотеза о нормальности отклонена, так как p-value < 0.05.\n",
    "- QQ-плот также показывает, что распределение далеко от нормального."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Применим логарифмирование для приближения распределения к нормальному"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tg = train['trip_duration'].apply(lambda x: np.log(x+1))\n",
    "sm.qqplot(tg)\n",
    "pv = shapiro(tg)[1]\n",
    "plt.title('p-value: ' + str(pv));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Логарифмирование не помогло. Распределение не стало нормальным"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Построение новых признаков\n",
    "### Объединим тестовые и обучающие выборки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.concat([train.drop(['dropoff_datetime','trip_duration'], axis=1),test], axis=0)\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Сгенерируем новые признаки из столбца pickup_datetime: \n",
    "    - hour - час посадки; \n",
    "    - dyear - день в году; \n",
    "    - dweek - день недели;\n",
    "    - weekofyear - неделя в году;\n",
    "    - minute - минута посадки;\n",
    "    - week_hour - день недели * 24 + час;\n",
    "    - dt - отклонение текущего времени посадки от минимального значения в секундах\n",
    "    - dweek_pass_count - день недели_число пассажиров. Разные дни могут характеризоваться разным кол-вом пассажиров.\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.pickup_datetime = pd.to_datetime(data.pickup_datetime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "data['hour'] = data.pickup_datetime.dt.hour\n",
    "data['dyear'] = data.pickup_datetime.dt.dayofyear\n",
    "data['dweek'] = data.pickup_datetime.dt.dayofweek\n",
    "data['weekofyear'] = data.pickup_datetime.dt.weekofyear\n",
    "data['minute'] = data['pickup_datetime'].dt.minute\n",
    "data['week_hour'] = data['dweek'] * 24 + data['hour']\n",
    "data['dt'] = (data.pickup_datetime - data.pickup_datetime.min()).dt.total_seconds()\n",
    "data['dweek_pass_count'] = list(map(lambda x,y: x*y, data['dweek'], data['passenger_count']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Кластеризуем координаты посадки и высадки\n",
    "- число кластеров 100 дало существенное улучшение качества модели"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_clust=100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#kmeans = KMeans(n_clusters=n_clust, random_state=0)\n",
    "kmeans = MiniBatchKMeans(n_clusters=n_clust, batch_size=100, random_state=0)\n",
    "kmeans.fit(data[['pickup_longitude','pickup_latitude']])\n",
    "cx = [c[0] for c in kmeans.cluster_centers_]\n",
    "cy = [c[1] for c in kmeans.cluster_centers_]\n",
    "data['pickup_cluster'] = kmeans.predict(data[['pickup_longitude','pickup_latitude']])\n",
    "data['dest_cluster'] = kmeans.predict(data[['dropoff_longitude', 'dropoff_latitude']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,8))\n",
    "plt.ylim(40.58, 40.93)\n",
    "plt.xlim(-74.1,-73.7)\n",
    "\n",
    "#colors = [np.random.rand(3,1) for i in range(n_clust)]\n",
    "cm = plt.get_cmap('gist_rainbow')\n",
    "colors = [cm(2.*i/n_clust) for i in range(n_clust)]\n",
    "colored = [colors[k] for k in data['pickup_cluster']]\n",
    "\n",
    "plt.scatter(data.pickup_longitude, data.pickup_latitude, s=0.01, color = colored)\n",
    "plt.scatter(cx, cy, s=20, color='black')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- При отображении координат посадки на графике и раскрашивании их по кластерам можно сделать следующие выводы:\n",
    "    - Наибольшее число посадок в такси приходит на остров Манхэттен, скорее всего из-за того, что исследуются поездки только  \"желтого такси\" NYC Yellow Cab, за которым закреплен район Манхэттен.\n",
    "    - В целом кластеры выделяются достаточно четко, что говорит об удачном выборе числа кластеров.\n",
    "    - Кластеризацию точек высадки пассажиров можно произвести по модели, обученной на точках посадки, так как большинство координат лежат находятся в одном и том же ареале."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Кодирование категориальных признаков\n",
    "- Кластеры координат посадки и высадки, а также признак день недели будем кодировать средним значением целевой переменной (trip_duration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target = train['trip_duration']\n",
    "train_df = pd.concat([data.iloc[:1444051],target], axis=1)\n",
    "train_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['pickup_cluster_mean'] = 0\n",
    "data['dest_cluster_mean'] = 0\n",
    "for cluster in tqdm.tqdm_notebook(range(n_clust)):\n",
    "    data.at[data['pickup_cluster']==cluster,'pickup_cluster_mean'] = \\\n",
    "    train_df[train_df['pickup_cluster']==cluster]['trip_duration'].mean()\n",
    "    \n",
    "    data.at[data['dest_cluster']==cluster,'dest_cluster_mean'] = \\\n",
    "    train_df[train_df['dest_cluster']==cluster]['trip_duration'].mean()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(8,6))\n",
    "ax = data['pickup_cluster_mean'].hist(bins=100, color='orange')\n",
    "ax = data['dest_cluster_mean'].hist(bins=100)\n",
    "ax.legend(['pickup_cluster_mean', 'dest_cluster_mean'], fontsize=15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Также закодируем средним категориальный признак dweek\n",
    "data['dweek_mean'] = 0\n",
    "for cluster in tqdm.tqdm_notebook(range(7)):\n",
    "    data.at[data['dweek']==cluster,'dweek_mean'] = \\\n",
    "    train_df[train_df['dweek']==cluster]['trip_duration'].mean()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "categor_features = [ 'pickup_cluster_mean', 'dest_cluster_mean', 'dweek_mean']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "- Координаты высадки имеют больше значений в правом хвосте, это говорит о том, что чаще посадка совершается в центре, и затем такси отвозит клеента в переферийные районы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Новые признаки после PCA преобразования координат. \n",
    "- PCA позволяет зеркально отразить координаты, что позволяет получить новые признаки из значений координат."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = np.vstack((data[['pickup_latitude', 'pickup_longitude']].values,\n",
    "                    data[['dropoff_latitude', 'dropoff_longitude']].values))\n",
    "\n",
    "pca = PCA().fit(coords)\n",
    "data['pickup_pca0'] = pca.transform(data[['pickup_latitude', 'pickup_longitude']])[:, 0]\n",
    "data['pickup_pca1'] = pca.transform(data[['pickup_latitude', 'pickup_longitude']])[:, 1]\n",
    "data['dropoff_pca0'] = pca.transform(data[['dropoff_latitude', 'dropoff_longitude']])[:, 0]\n",
    "data['dropoff_pca1'] = pca.transform(data[['dropoff_latitude', 'dropoff_longitude']])[:, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Признаки расстояний между точками высадки и посадки\n",
    "- Расстояния между точками на сфере можно рассчитать несколькими способами:\n",
    "    - Формула гаверсинусов (haversine) для небольших расстояний на сфере.\n",
    "    - Сферическая теорема косинусов.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def haversine_array(lat1, lng1, lat2, lng2):\n",
    "    lat1, lng1, lat2, lng2 = map(np.radians, (lat1, lng1, lat2, lng2))\n",
    "    AVG_EARTH_RADIUS = 6371  # in km\n",
    "    lat = lat2 - lat1\n",
    "    lng = lng2 - lng1\n",
    "    d = np.sin(lat * 0.5) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(lng * 0.5) ** 2\n",
    "    h = 2 * AVG_EARTH_RADIUS * np.arcsin(np.sqrt(d))\n",
    "    return h\n",
    "\n",
    "def cos_dist(lat1, lng1, lat2, lng2):\n",
    "    lat1, lng1, lat2, lng2 = map(np.radians, (lat1, lng1, lat2, lng2))\n",
    "    AVG_EARTH_RADIUS = 6371  # in km\n",
    "    lat = lat2 - lat1\n",
    "    lng = lng2 - lng1\n",
    "    d = np.sin(lat1)*np.sin(lat2) + np.cos(lng1) * np.cos(lng2) * np.cos(lng)\n",
    "    h = AVG_EARTH_RADIUS * np.arccos(d)\n",
    "    return h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['distance_haversine'] = haversine_array(data['pickup_latitude'].values, data['pickup_longitude'].values, data['dropoff_latitude'].values, data['dropoff_longitude'].values)\n",
    "data['cos_dist'] = cos_dist(data['pickup_latitude'].values, data['pickup_longitude'].values, data['dropoff_latitude'].values, data['dropoff_longitude'].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Манхеттенские расстояния для обычных координат и отраженных координат с помощью PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['manhattan'] = np.abs(data['dropoff_latitude'] - data['pickup_latitude']) + np.abs(data['dropoff_longitude'] - data['pickup_longitude'])\n",
    "data['pca_manhattan'] = np.abs(data['dropoff_pca1'] - data['pickup_pca1']) + np.abs(data['dropoff_pca0'] - data['pickup_pca0'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Дополнительные признаки - центр на прямой линии между точками посадки и высадки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['center_latitude'] = (data['pickup_latitude'].values + data['dropoff_latitude'].values) / 2\n",
    "data['center_longitude'] = (data['pickup_longitude'].values + data['dropoff_longitude'].values) / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Признаки: погода\n",
    "- Информация о [погоде в Центральном парке за 2016 год](https://www.kaggle.com/mathijs/weather-data-in-new-york-city-2016/data) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_data = pd.read_csv('../../data/weather_data_nyc_centralpark_2016(1).csv')\n",
    "w_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_data.date = pd.to_datetime(w_data.date)\n",
    "# Преобразуем значение T  в признаках - snow fall и precipitation\n",
    "falls = [ 0.01 if c=='T' else float(c) for c in w_data['snow fall']]\n",
    "rain = [ 0.01 if c=='T' else float(c) for c in w_data['precipitation']]\n",
    "w_data['snow fall']= falls\n",
    "w_data['precipitation'] = rain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_data['dyear'] = w_data.date.dt.dayofyear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_data['snow depth'] = w_data['snow depth'].apply(lambda x: 1 if x=='T' else float(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in w_data.columns[3:-1]:\n",
    "    data[col] = 0.0\n",
    "    print(col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for day in tqdm.tqdm_notebook(data['dyear'].value_counts().keys(), total = data['dyear'].value_counts().shape[0]):\n",
    "    for col in w_data.columns[3:-1]:\n",
    "        data.at[data['dyear'] == day, col] = w_data[w_data['dyear'] == day][col].values[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Для обучения будем использовать следующие признаки погоды:\n",
    "- average temperature - средняя температура за день,\n",
    "- precipitation - количество осадков в виде дождя,\n",
    "- snow fall\t - количество осадков в виде снега,\n",
    "- snow depth - глубина снежного покрова"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Дискретизируем признак dyear и добавим его как новый в основную выборку "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "data['dyear'].hist(bins=200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- На распределении признака  dyear видно 19 четких кластеров.\n",
    "- Выделим эти кластеры и закодируем средним по целевой переменной."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_f = list(map(lambda x: list(map(int, re.findall(r'\\d+', x))), pd.qcut(data.dyear, 19).values.astype(str)))\n",
    "data['new_dyear'] = list(map(lambda x: mean([x[0], x[2]]), new_f))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "data['new_dyear'].hist(bins=200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Также закодируем средним категориальный признак dweek\n",
    "target = train['trip_duration']\n",
    "train_df = pd.concat([data.iloc[:1444051],target], axis=1)\n",
    "\n",
    "data['new_dyear_mean'] = 0\n",
    "for cluster in tqdm.tqdm_notebook(data['new_dyear'].value_counts().keys()):\n",
    "    data.at[data['new_dyear']==cluster,'new_dyear_mean'] = \\\n",
    "    train_df[train_df['new_dyear']==cluster]['trip_duration'].mean()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Преобразование бинарных признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.store_and_fwd_flag = data.store_and_fwd_flag.apply(lambda x: 1 if x=='N' else 0)\n",
    "data.store_and_fwd_flag.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['vendor_id'] = data['vendor_id'].apply(lambda x: x-1)\n",
    "data['vendor_id'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='blue'>Сохраним данные</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.id = data.id.apply(lambda x: int(re.findall(r'\\d+', x)[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.to_hdf('data.h5','df')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train['trip_duration'].to_hdf('target.h5','df')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='green'>Считаем данные</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target = pd.read_hdf('target.h5','df')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_hdf('data.h5','df')\n",
    "\n",
    "data.shape, target.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  4. Визуальный анализ всех признаков \n",
    "#### Анализ распределений признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "real_features = ['passenger_count',\n",
    "       'pickup_longitude', 'pickup_latitude', 'dropoff_longitude',\n",
    "       'dropoff_latitude', 'hour', 'dyear',\n",
    "       'weekofyear', 'minute', 'week_hour', 'dt', 'dweek_pass_count',\n",
    "                'pickup_pca0', 'pickup_pca1',\n",
    "       'dropoff_pca0', 'dropoff_pca1', 'distance_haversine', 'cos_dist',\n",
    "       'pca_manhattan', 'manhattan', 'center_latitude', 'center_longitude',\n",
    "       'average temperature', 'precipitation', 'snow fall', 'snow depth']\n",
    "categor_features = ['pickup_cluster_mean', 'dest_cluster_mean', 'dweek_mean', 'new_dyear_mean']\n",
    "bin_features = ['vendor_id','store_and_fwd_flag']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_names = real_features + categor_features + bin_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(feature_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sns.set_style(\"darkgrid\")\n",
    "fig, axes = plt.subplots(nrows=6, ncols=6, figsize=(20,20))\n",
    "fig.subplots_adjust(hspace=0.5) # Увеличиваем пространство между subplots\n",
    "for idx, f in enumerate(feature_names):\n",
    "    ax = axes[int(idx/6),idx%6]\n",
    "    data[f].hist(ax=ax, bins=100)\n",
    "    #ax.legend([f], fontsize=10, loc='upper left')\n",
    "    ax.set_title(f, fontsize=15)\n",
    "    #ax.set_xlabel('index')\n",
    "    #ax.set_ylabel(u'Число сессий')\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Судя по распределениям координатные признаки и признаки на их основе имеют большие выбросы.\n",
    "- Выбросы по координатам объясняются удаленными районами города, куда таксистам приходится ездить. В то же время основная масса поездок приходится на центральные районы (Манхеттен, Бруклин, аэропорт). \n",
    "- Аналогичная ситуация наблюдается с другими признаками, являющимися производными от координатных признаков.\n",
    "- Координаты высадки 'dest_cluster_mean' имеют больше значений в правом хвосте по сравнению с координатами посадки 'pickup_cluster_mean', это говорит о том, что чаще посадка совершается в центре, и затем такси отвозит клеента в переферийные районы.\n",
    "- Можно проверить фильтрацию координат посадки и высадки - заменить координаты удаленных районов медианными значениями.\n",
    "- Кластеризация по координатам позволяет улучшить распределение, сделать его менее скошенным и вариативным. Для линейной модели это может быть критично, но для случайного леса нет.\n",
    "- Бинарный признак 'store_and_fwd_flag' (флаг задержки в передачи данных о поездке вендору) очень несбалансированный, практически константный. В большинстве случаев задержек в передачи данных не было. Скорее всего он не будет вносить заметный вклад в улучшение качества прогноза.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Корреляционный анализ признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_df = data.iloc[:target.shape[0]][feature_names]\n",
    "corr_df['trip_duration'] = target.values\n",
    "corr = corr_df[feature_names+['trip_duration']].corr(method='pearson')\n",
    "# Генерируем маску для отсечения верхней диагонали матрицы\n",
    "mask = np.zeros_like(corr, dtype=np.bool)\n",
    "mask[np.triu_indices_from(mask)] = True # Верхнему треугольнику присваиваем True\n",
    "\n",
    "plt.figure(figsize=(16,16))\n",
    "sns.heatmap(corr, mask=mask, vmax=1, center=0, annot=True, fmt='.1f',\n",
    "            square=True, linewidths=.5, cbar_kws={\"shrink\": 1});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- На матрице корреляций по Пирсону видны пары признаки с высокой, единичной корреляцией. Что говорит о линейной зависисимости этих признаков. \n",
    "- Если избавиться по одному из признаков из каждой такой пары, то можно избежать переобучение на линейных моделях. Но скорее всего линейные модели в этой задаче будут работать гораздо хуже моделей случайного леса и крадиентного бустинга из-за того, что признаки имеют плохие несимметричные распределения с большим количеством выбросов. Но если бы линейная модель была основной, то можно было избавиться от следующих признаков: **'center_longitude', 'pca_manhattan', 'pickup_pca1','dropoff_pca1', 'week_hour', 'dyear'**.\n",
    "- Также из корреляционной матрицы видно, что нет признаков, которые имеют единичную положительную корреляцию с целевой переменной trip_duration. Наличие такого признака делало бы бессмысленным всю работу по проектированию предсказательной модели.\n",
    "- В то же время видно, что функции расстояния - гаверсинус ('distance_haversine') и манхеттенская метрика ('pca_manhattan') достаточно сильно коррелируют с целевой переменной, что вполне естественно, поскольку целевая переменная  - это продолжительность поездки в секундах.\n",
    "- Расстояние расчитанное по сферической теореме косинусов ('cos_dist') не коррелирует с целевой переменной, что может говорить о недостатке этой мере расстояния. [Из википедии:](https://en.wikipedia.org/wiki/Great-circle_distance)В случае маленьких расстояний и небольшой разрядности вычисления (количество знаков после запятой), использование формулы может приводить к значительным ошибкам связанным с округлением. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Выбор метрики\n",
    "В качестве метрики выбрана средняя абсолютная ошибка (MAE):\n",
    "    - MAE - наиболее удобная метрика для задач регрессии, показывает отклонение в единицах измерения целевой переменной."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Выбор модели\n",
    "Рассмотрим несколько типов моделей:\n",
    "- Случайный лес - подходит в нашем случае, поскольку восстанавливаемая функция скорее всего нелинейная. Также распределения многих признаков далеки от нормальных со множеством выбросов. С такими нетипичными распределениями лучше справляются модели основанные на ансамблях алгоритмов.\n",
    "- Линейную модель LASSO (L1-регуляризация) попробуем применить для отбора признаков.\n",
    "- Градиентный бустинг LightGBM - быстрая версия бустинга от Microsoft. Применение быстрой модели позволит на слабой машине проверить Стекинг - способ ансамблирования моделей с построением метапризнаков. Случайный лес со стекингом для конкретной задачи (~1,5 млн записей) в несколько раз увеличивает время расчета."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Типы признаков\n",
    "- В обучении модели не будем использовать признаки: 'dyear', 'dweek', 'pickup_cluster' и 'dest_cluster', так как мы их заменили закодированными средним значением целевой переменной.\n",
    "- Основную выборку будем делить на обучающую и валидационную в соотношении 70/30 с перемешиванием, так как мы имеем почти 1.5 млн. записей, следовательно, 70% вполне хватит на обучение модели, а 30% - хороший объем для валидации."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "real_features = ['passenger_count',\n",
    "       'pickup_longitude', 'pickup_latitude', 'dropoff_longitude',\n",
    "       'dropoff_latitude', 'hour', \n",
    "       'weekofyear', 'minute', 'week_hour', 'dt', 'dweek_pass_count',\n",
    "                'pickup_pca0', 'pickup_pca1',\n",
    "       'dropoff_pca0', 'dropoff_pca1', 'distance_haversine', 'cos_dist',\n",
    "       'pca_manhattan', 'manhattan', 'center_latitude', 'center_longitude',\n",
    "       'average temperature', 'precipitation', 'snow fall', 'snow depth']\n",
    "categor_features = ['pickup_cluster_mean', 'dest_cluster_mean', 'dweek_mean', 'new_dyear_mean']\n",
    "bin_features = ['vendor_id','store_and_fwd_flag']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RandomForestRegressor\n",
    "- Возьмем в качестве базовой модели Случайный лес"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_names = real_features + categor_features + bin_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_full = data[feature_names].values\n",
    "X=X_full[:1444051]\n",
    "X_tst = X_full[1444051:]\n",
    "y = target.values\n",
    "X.shape, y.shape, X_tst.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Разобъем выборку на обучающую и валидационную\n",
    "X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.3, \n",
    "                                                     random_state=RANDOM_STATE)\n",
    "X_train.shape, X_valid.shape, y_train.shape, y_valid.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kfold = KFold(n_splits = 5, shuffle = True, random_state = RANDOM_STATE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "rf = RandomForestRegressor(n_estimators = 30,\n",
    "                           max_depth = 10,\n",
    "                           random_state=RANDOM_STATE)\n",
    "mean_score = cross_val_score(rf, X_train, y_train, cv=kfold, n_jobs=-1, scoring = 'neg_mean_absolute_error').mean()\n",
    "print ('mae: {}'.format(mean_score))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Случайный лес работает слишком долго на таком объеме данных.\n",
    "- Проверим с помощью линейной модели LASSO значимость признаков и затем попробуем применить отобранные значимые признаки на базовой модели случайного леса."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Масштабирование признаков"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Вещественные признаки\n",
    "scaler = StandardScaler()\n",
    "X_full = scaler.fit_transform(data[real_features+categor_features])\n",
    "X_full.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Бинарные признаки\n",
    "X_full = np.concatenate((X_full,data[bin_features].values), axis=1)\n",
    "X_full.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X=X_full[:1444051]\n",
    "X_tst = X_full[1444051:]\n",
    "y = target.values\n",
    "X.shape, y.shape, X_tst.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_scaled_train, X_scaled_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=RANDOM_STATE)\n",
    "X_scaled_train.shape, X_scaled_test.shape, y_train.shape, y_test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Линейная модель\n",
    "#### Анализ значимости признаков с помощью модели Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kfold = KFold(n_splits = 5, shuffle = True, random_state = RANDOM_STATE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "lasso = LassoCV(cv = kfold, n_jobs=-1)\n",
    "lasso.fit(X_scaled_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_absolute_error(y_test, lasso.predict(X_scaled_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_names = real_features + categor_features + bin_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Функция, возвращающая веса линейной модели, убывающие по абсолютному значению\n",
    "def getCoefficients(model):    \n",
    "    coefs = pd.DataFrame(model.coef_, feature_names)\n",
    "    coefs.columns = [\"coef\"]\n",
    "    coefs[\"abs\"] = coefs.coef.apply(np.abs)\n",
    "    return coefs.sort_values(by=\"abs\", ascending=False).drop([\"abs\"], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "getCoefficients(lasso)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "good_features = ['distance_haversine', 'dest_cluster_mean', 'pickup_cluster_mean', 'dweek_mean', \n",
    "                 'new_dyear_mean', 'hour', 'center_latitude', 'manhattan']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Проверим как изменилась ошибка базового случайного леса на отобранных значимых признаках good_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_full = data[good_features].values\n",
    "X=X_full[:1444051]\n",
    "X_tst = X_full[1444051:]\n",
    "y = target.values\n",
    "# Разобъем выборку на обучающую и валидационную\n",
    "X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.3, \n",
    "                                                     random_state=RANDOM_STATE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "rf = RandomForestRegressor(n_estimators = 30,\n",
    "                           max_depth = 10,\n",
    "                           random_state=RANDOM_STATE)\n",
    "mean_score = cross_val_score(rf, X_train, y_train, cv=kfold, n_jobs=-1, scoring = 'neg_mean_absolute_error').mean()\n",
    "print ('mae: {}'.format(mean_score))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Lasso метод выкинул значащие признаки, так как MAE повысилась после обучения на неполной выборке. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Попробуем улучшить качество за счет фильтрация выбросов некоторых признаков\n",
    "- Заменим едианными значениями выбросы по координатам"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "plt1 = plt.scatter(range(data.shape[0]), np.sort(data['pickup_longitude'].apply(np.abs)))\n",
    "plt2 = plt.scatter(range(data.shape[0]), np.sort(data['dropoff_longitude'].apply(np.abs)), color='orange')\n",
    "plt.xlabel('index')\n",
    "plt.legend([plt1, plt2], ['pickup_longitude', 'dropoff_longitude'], fontsize=15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.at[data['pickup_longitude'] < -85.0, 'pickup_longitude'] = data['pickup_longitude'].median()\n",
    "data.at[data['dropoff_longitude'] < -85.0, 'dropoff_longitude'] = data['pickup_longitude'].median()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "plt1 = plt.scatter(range(data.shape[0]), np.sort(data['pickup_latitude']))\n",
    "plt2 = plt.scatter(range(data.shape[0]), np.sort(data['dropoff_latitude']), color='orange') \n",
    "plt.xlabel('index')\n",
    "plt.legend([plt1, plt2], ['pickup_latitude', 'dropoff_latitude'], fontsize=15);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.at[data['pickup_latitude'] > 45.0, 'pickup_latitude'] = data['pickup_latitude'].median()\n",
    "data.at[data['dropoff_latitude'] > 45.0, 'dropoff_latitude'] = data['dropoff_latitude'].median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Обучим базовую модель RandomForest на отфильтрованных данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_names = real_features + categor_features + bin_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_full = data[feature_names].values\n",
    "X=X_full[:1444051]\n",
    "X_tst = X_full[1444051:]\n",
    "y = target.values\n",
    "# Разобъем выборку на обучающую и валидационную\n",
    "X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.3, \n",
    "                                                     random_state=RANDOM_STATE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "rf = RandomForestRegressor(n_estimators = 30,\n",
    "                           max_depth = 10,\n",
    "                           random_state=RANDOM_STATE)\n",
    "mean_score = cross_val_score(rf, X_train, y_train, cv=kfold, n_jobs=-1, scoring = 'neg_mean_absolute_error').mean()\n",
    "print ('mae: {}'.format(mean_score))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Ошибка практически не изменилась по сравнению с результатом на неотфильтрованных данных (mae: -202.35252728664412)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LightGBM\n",
    "- Попробуем улучшить результат за счет применения градиентного бустинга\n",
    "- Будем использовать реализацию бустинга от Microsoft - LightGBM, которая значительно быстрее традиционного XGBoost. Использование быстрой модели позволит обучать множество вариантов моделей с различными гиперпапраметрами и усреднять результаты."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "real_features = ['passenger_count',\n",
    "       'pickup_longitude', 'pickup_latitude', 'dropoff_longitude',\n",
    "       'dropoff_latitude', 'hour', #'dyear', #'dweek',\n",
    "       'weekofyear', 'minute', 'week_hour', 'dt', 'dweek_pass_count',\n",
    "                'pickup_pca0', 'pickup_pca1',\n",
    "       'dropoff_pca0', 'dropoff_pca1', 'distance_haversine', 'cos_dist',\n",
    "       'manhattan', 'center_latitude', 'center_longitude', 'pca_manhattan', \n",
    "       'average temperature', 'precipitation', 'snow fall', 'snow depth'\n",
    "                ]\n",
    "categor_features = ['pickup_cluster_mean', 'dest_cluster_mean', 'dweek_mean', 'new_dyear_mean']\n",
    "bin_features = ['vendor_id','store_and_fwd_flag']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Деление выборки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_names = real_features + categor_features + bin_features\n",
    "#feature_names = good_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_full = data[feature_names].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X=X_full[:1444051]\n",
    "X_tst = X_full[1444051:]\n",
    "y = target.values\n",
    "X.shape, y.shape, X_tst.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.3, \n",
    "                                                     random_state=7)\n",
    "X_train.shape, X_valid.shape, y_train.shape, y_valid.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Модель LightGBM с праметрами по умолчанию"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = lgb.LGBMRegressor(\n",
    "            seed=RANDOM_STATE,\n",
    "            objective='regression',\n",
    "            nthread=-1\n",
    "        )\n",
    "model.fit(X_train, y_train)        \n",
    "mae = metrics.mean_absolute_error(y_valid, model.predict(X_valid))\n",
    "print('MAE:{}'.format(mae))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Пока качество хуже случайного леса \"из коробки\".\n",
    "- Но попробуем настроить гиперпараметры."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Настройка гиперпараметров модели с помошью библиотеки hyperopt\n",
    "- Выборка очень большая, более 1 млн. строк, поэтому обойдемся без кросс-валидации. Валидироваться будем просто на отложенной выборке один раз."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hyperopt import hp\n",
    "from hyperopt import fmin, Trials, tpe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hyperopt_objective(params):\n",
    "    model = lgb.LGBMRegressor(\n",
    "            seed=RANDOM_STATE,\n",
    "            n_estimators = params['num_est'],\n",
    "            max_depth = params['max_dpth'],\n",
    "            num_leaves = params['num_lvs'],\n",
    "            learning_rate = params['lrn_rate'],\n",
    "            reg_alpha = params['reg_alph'],\n",
    "            reg_lambda = params['reg_lmbd'],\n",
    "            objective='regression',\n",
    "            nthread=-1\n",
    "        )    \n",
    "    model.fit(X_train, y_train)\n",
    "    mae = metrics.mean_absolute_error(y_valid, model.predict(X_valid))\n",
    "    return mae"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_space = {\n",
    "    'num_est': hp.choice('num_est', (4000,4500,5000,5500,6000)), \n",
    "    'max_dpth': hp.choice('max_dpth', (30,33,35,38)),\n",
    "    'num_lvs': hp.choice('num_lvs', (15,18,20,23,27)),\n",
    "    'lrn_rate': hp.uniform('lrn_rate', 0.2, 0.4),\n",
    "    'reg_alph': hp.uniform('reg_alph', 10, 50),\n",
    "    'reg_lmbd': hp.uniform('reg_lmbd', 10, 50)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "trials = Trials()\n",
    "best_params = fmin(\n",
    "    hyperopt_objective,\n",
    "    space=params_space,\n",
    "    algo=tpe.suggest,\n",
    "    max_evals=100,\n",
    "    trials=trials,\n",
    "    verbose = 2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Сохраним лог hyperopt\n",
    "with open('hyperopt_trials.pkl', 'wb') as hyp_trl:\n",
    "    pickle.dump(trials.trials, hyp_trl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Пример лога поиска оптимальных параметров с помощью hyperopt:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trials.trials[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### График ошибки на отложенной выборке, который показывает как hyperopt подбирал параметры"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "plt.plot(range(len(trials.losses())), trials.losses());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "- Судя по графику поиск происходил случайным образом"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min(trials.losses())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Обучение на оптимальных гиперпараметрах"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = lgb.LGBMRegressor(\n",
    "            seed=RANDOM_STATE,\n",
    "            n_estimators=5500,\n",
    "            max_depth=38,\n",
    "            num_leaves=27,\n",
    "            learning_rate=0.2029823209429007,\n",
    "            reg_alpha=43.63743528388524,\n",
    "            reg_lambda=42.56637567284501,\n",
    "            objective='regression',\n",
    "            nthread=-1\n",
    "        )\n",
    "model.fit(X_train, y_train)        \n",
    "mae = metrics.mean_absolute_error(y_valid, model.predict(X_valid))\n",
    "print('MAE:{}'.format(mae))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Видим заметное улучшение качества прогнозирования за счет подбора потимальных гиперпараметров модели с помощью бибилиотеки hyperopt по сравнению с дeфолтной моделью lightgbm (MAE:265.562706614055)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 9. Построение валидационных, тестовых и тренировочных кривых.\n",
    "- Перебор 100 комбинаций гиперпараметров lightgbm занял около 8 часов. Столь длительное время объясняется главным образом количеством деревьев (от 4000 до 6000). Построение же кросс-валидационных кривых для такого количества деревьев может занять в несколько раз больше времени. Поэтому было принято решение построить валидационные кривые только для параметров 'max_dpth' и 'num_lvs' для 1000 деревьев, и оценить по ним как обучается модель."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kfold = KFold(n_splits = 3, shuffle = True, random_state = RANDOM_STATE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Кривые валидации, тестовая и тренировочная для разных значений параметра ***max_depth***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_mae = []\n",
    "valid_mae = []\n",
    "test_mae = [] \n",
    "max_depth_grid = [10,40,60,80,100]\n",
    "\n",
    "# Обучаем на тренировочном датасете\n",
    "for depth in tqdm.tqdm_notebook(max_depth_grid):\n",
    "    model = lgb.LGBMRegressor(\n",
    "            seed=RANDOM_STATE,\n",
    "            n_estimators=1000,\n",
    "            max_depth=depth,\n",
    "            num_leaves=27,\n",
    "            learning_rate=0.2029823209429007,\n",
    "            reg_alpha=43.63743528388524,\n",
    "            reg_lambda=42.56637567284501,\n",
    "            objective='regression',\n",
    "            nthread=-1\n",
    "        )\n",
    "    temp_train_mae = []\n",
    "    temp_valid_mae = []\n",
    "    temp_test_mae = []\n",
    "    for train_index, valid_index in kfold.split(X_train, y_train):\n",
    "        X_trn, X_vld = X_train[train_index], X_train[valid_index]\n",
    "        y_trn, y_vld = y_train[train_index], y_train[valid_index]\n",
    "        model.fit(X_trn, y_trn)\n",
    "        temp_train_mae.append(metrics.mean_absolute_error(y_trn, model.predict(X_trn)))\n",
    "        temp_valid_mae.append(metrics.mean_absolute_error(y_vld, model.predict(X_vld)))\n",
    "    train_mae.append(temp_train_mae)\n",
    "    valid_mae.append(temp_valid_mae)\n",
    "    test_mae.append(metrics.mean_absolute_error(y_valid, model.predict(X_valid)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_cv_trn_tst(train_mae, valid_mae, test_mae, grid, label):\n",
    "    train_mae, valid_mae, test_mae = np.asarray(train_mae), np.asarray(valid_mae), np.asarray(test_mae)\n",
    "    plt.style.use('ggplot')\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(12, 8))\n",
    "    ax.plot(grid, train_mae.mean(axis=1), alpha=0.5, color='blue', label='train')\n",
    "    ax.plot(grid, valid_mae.mean(axis=1), alpha=0.5, color='red', label='cv')\n",
    "    ax.fill_between(grid, valid_mae.mean(axis=1) - valid_mae.std(axis=1), valid_mae.mean(axis=1) + valid_mae.std(axis=1), color='#888888', alpha=0.4)\n",
    "    ax.fill_between(grid, valid_mae.mean(axis=1) - 2*valid_mae.std(axis=1), valid_mae.mean(axis=1) + 2*valid_mae.std(axis=1), color='#888888', alpha=0.2)\n",
    "    ax.plot(grid, test_mae, alpha=0.5, color='orange', label='test')\n",
    "    ax.legend(loc='best', fontsize=15)\n",
    "    #ax.set_ylim([0.88,1.02])\n",
    "    ax.set_ylabel(\"MAE\", fontsize=15)\n",
    "    ax.set_xlabel(label, fontsize=15)\n",
    "    ax.yaxis.set_tick_params(labelsize=15)\n",
    "    ax.xaxis.set_tick_params(labelsize=15)\n",
    "    ax.set_title('CV, train, test-curves', fontsize=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_cv_trn_tst(train_mae, valid_mae, test_mae, max_depth_grid, 'max_depth')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Видно, что после  'max_depth' = 40 качество модели не улучшается. То есть, параметр, подобранный функцийе hyperopt является оптимальным."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Кривые валидации, тестовая и тренировочная для разных значений параметра ***num_leaves***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_mae = []\n",
    "valid_mae = []\n",
    "test_mae = [] \n",
    "num_lvs_grid = [40,50,55,60,65]\n",
    "\n",
    "# Обучаем на тренировочном датасете\n",
    "for lvs in tqdm.tqdm_notebook(num_lvs_grid):\n",
    "    model = lgb.LGBMRegressor(\n",
    "            seed=RANDOM_STATE,\n",
    "            n_estimators=1000,\n",
    "            max_depth=40,\n",
    "            num_leaves=lvs,\n",
    "            learning_rate=0.2029823209429007,\n",
    "            reg_alpha=43.63743528388524,\n",
    "            reg_lambda=42.56637567284501,\n",
    "            objective='regression',\n",
    "            nthread=-1\n",
    "        )\n",
    "    temp_train_mae = []\n",
    "    temp_valid_mae = []\n",
    "    temp_test_mae = []\n",
    "    for train_index, valid_index in kfold.split(X_train, y_train):\n",
    "        X_trn, X_vld = X_train[train_index], X_train[valid_index]\n",
    "        y_trn, y_vld = y_train[train_index], y_train[valid_index]\n",
    "        model.fit(X_trn, y_trn)\n",
    "        temp_train_mae.append(metrics.mean_absolute_error(y_trn, model.predict(X_trn)))\n",
    "        temp_valid_mae.append(metrics.mean_absolute_error(y_vld, model.predict(X_vld)))\n",
    "    train_mae.append(temp_train_mae)\n",
    "    valid_mae.append(temp_valid_mae)\n",
    "    test_mae.append(metrics.mean_absolute_error(y_valid, model.predict(X_valid)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_cv_trn_tst(train_mae, valid_mae, test_mae, num_lvs_grid, 'num_leaves')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Видно, что после значения параметра 'num_leaves' = 50 тестовая кривая перестает снижаться. Поэтому на этом значении 'num_leaves' стоит остановиться. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Проверим качество модели на отложенной выборке с количеством деревьев, найденных с помощью hyperopt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = lgb.LGBMRegressor(\n",
    "            seed=RANDOM_STATE,\n",
    "            n_estimators=5500,\n",
    "            max_depth=38,\n",
    "            num_leaves=50,\n",
    "            learning_rate=0.2029823209429007,\n",
    "            reg_alpha=43.63743528388524,\n",
    "            reg_lambda=42.56637567284501,\n",
    "            objective='regression',\n",
    "            nthread=-1\n",
    "        )\n",
    "model.fit(X_train, y_train)        \n",
    "mae = metrics.mean_absolute_error(y_valid, model.predict(X_valid))\n",
    "print('MAE:{}'.format(mae))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Качество немного улучшилось без риска переобучения.\n",
    "- Оставим эти гиперпараметры и попробуем применить Стекинг."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10. Стекинг\n",
    "- Идея стекинга заключается в генерировании метапризнаков с помощью некоторой базовой модели и использование полученных метапризнаков совместно с основной выборкой в итоговой метамодели: [Методы ансамблирования обучающихся алгоритмов](http://www.machinelearning.ru/wiki/images/5/56/Guschin2015Stacking.pdf)\n",
    "- В текущем проекте использовался модифицированный вариант стекинга, когда деление исходной выборки осуществлялось методом bootstrap.\n",
    "- В качестве метамодели и основной модели использовалась LightGBM.\n",
    "- Из-за большого объема исходных данных стекинг на большом количестве bootstrap-выборок занимает очень много времени на среднем по производительности PC. В то же время даже на трех bootstrap-выборках удалось добиться существенного снижения ошибки. \n",
    "- Во время обучения стекинг осуществлялся в 2 прохода для разных значений random_state, а полученные результаты 2-х моделей усреднялись. Такой подход также привел к снижению ошибки."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Разбиение с помощью bootstrap\n",
    "def get_bootstrap_samples(X, y, n_samples):  \n",
    "    X_train_1 = []\n",
    "    y_train_1 = []\n",
    "    X_train_2 = []\n",
    "    y_train_2 = []\n",
    "    # генерируем матрицу индексов подвыборок из генеральной выборки\n",
    "    indices = np.random.randint(0, y.size, (n_samples, y.size)) \n",
    "    # массив выборок\n",
    "    j=0\n",
    "    for i in range(int(n_samples/2)):\n",
    "        X_train_1.append(X[indices[j]])\n",
    "        y_train_1.append(y[indices[j]])\n",
    "        X_train_2.append(X[indices[j+1]])\n",
    "        y_train_2.append(y[indices[j+1]])  \n",
    "        j += 2\n",
    "    return X_train_1, y_train_1, X_train_2, y_train_2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Модель LightGBM с предварительно настроенными гиперпараметрами\n",
    "def lgbm_model():\n",
    "    model = lgb.LGBMRegressor(\n",
    "            seed=RANDOM_STATE,\n",
    "            n_estimators=5500,\n",
    "            max_depth=38,\n",
    "            num_leaves=27,\n",
    "            learning_rate=0.2029823209429007,\n",
    "            reg_alpha=43.63743528388524,\n",
    "            reg_lambda=42.56637567284501,\n",
    "            objective='regression',\n",
    "            nthread=-1\n",
    "        )\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "KK = [3] # количество выборок bootstrap\n",
    "y_pred_best = None\n",
    "meta_mae_0 = inf\n",
    "\n",
    "basic_mae = []\n",
    "meta_mae = []\n",
    "y_pred = []\n",
    "seed = [42,7] # Варианты сидов\n",
    "n_samples = len(seed)\n",
    "for i in tqdm.tqdm_notebook(range(n_samples)):\n",
    "    np.random.seed(seed[i])\n",
    "    kk = KK[0] \n",
    "    X_train_1, y_train_1, X_train_2, y_train_2 = get_bootstrap_samples(X, y, kk*2)\n",
    "    \n",
    "    for k in tqdm.tqdm_notebook(range(kk)):\n",
    "        # Обучим базовую модель на выборке X_train_1, y_train_1\n",
    "        bst = lgbm_model()\n",
    "        bst.fit(X_train_1[k], y_train_1[k])\n",
    "                \n",
    "        # make prediction\n",
    "        basic_mae.append(metrics.mean_absolute_error(y_valid, bst.predict(X_valid)))   \n",
    "\n",
    "        # Предсказание обученного алгоритма для выборок X_train_2 и X_valid. Метапризнаки\n",
    "        MF2 = bst.predict(X_train_2[k])\n",
    "        MF0 = bst.predict(X_valid)\n",
    "        MF_tst = bst.predict(X_tst)\n",
    "\n",
    "        # Получим новую обучающую выборку (объединим X_train_2 c метапризнаком MF2)\n",
    "        new_train = np.concatenate((X_train_2[k],MF2.reshape(MF2.size, 1)), axis=1)\n",
    "\n",
    "        # Получим новую валидационную выборку (объединим X_valid c метапризнаком MF0)\n",
    "        new_valid = np.concatenate((X_valid,MF0.reshape(MF0.size, 1)), axis=1)\n",
    "\n",
    "        # Получим новую тестовую выборку (объединим X_tst c метапризнаком MF_tst)\n",
    "        new_tst = np.concatenate((X_tst,MF_tst.reshape(MF_tst.size, 1)), axis=1)\n",
    "\n",
    "        # Обучим метамодель        \n",
    "        meta_model = lgbm_model()\n",
    "        meta_model.fit(new_train, y_train_2[k])\n",
    "\n",
    "        meta_mae.append(metrics.mean_absolute_error(y_valid, meta_model.predict(new_valid)))\n",
    "        \n",
    "        y_pred.append(meta_model.predict(new_tst))\n",
    "\n",
    "basic_mean = mean(basic_mae)\n",
    "meta_mean = mean(meta_mae)\n",
    "\n",
    "print ('basic_mean: {}, meta_mean: {}'.format(basic_mean, meta_mean))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('{}%'.format(round((100*(155 - meta_mean)/155),0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Стекинг позволил снизить MAE на 20% (123.9) по сравнению с моделью LightGBM с просто настроенными параметрами."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Подготовка данных для отправки на kaggle.com"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = np.vstack(y_pred).T\n",
    "predictions = predictions.mean(axis=1)\n",
    "test['trip_duration'] = predictions.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Проверка на NAN\n",
    "test['trip_duration'].isnull().values.any()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test['trip_duration'] = test['trip_duration'].apply(lambda x: x*(-1) if x<0 else x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load to csv\n",
    "test[['id','trip_duration']].to_csv('submit_.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Резултат на kaggle (метрика Root Mean Squared Logarithmic Error):\n",
    "- Private score: 0.39744 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11. Выводы.\n",
    "- В процессе анализа признаков было сгенерировано большое количество временных и гео-признаков. Анализ этих признаков позволяет получить много информации об особенностях транспортных пассажирских перевозках в исследуемой локации. То есть, помимо решенной задачи прогнозирования поездки, эта работа позволяет сделать выводы по дальнейшему улучшению качества поездок.\n",
    "- Далеко не все признаки сделали существенный вклад в итоговый прогноз. Но это не всегда говорит о том, что признак бесполезен. Часто грамотное кодирование признака и применение подходящей модели позволяет извлечь из него много полезной информации.\n",
    "- Настройка гиперапраметров модели LightGBM позволило получить результат существенно лучше случайного леса и линейных моделей.\n",
    "- Применение технологии стекинга (обучение на метапризнаках) позволило снизить ошибку на 20% по сравнению с моделью, для которой осуществлялась только настройка параметров.\n",
    "- Дальнейшее развитие модели:\n",
    "    - Так как в задаче используются распределенные в пространстве данные, то можно представить признаковое пространство в виде геораспределенной сетки и применить сверточную сеть по аналогии с изображением.\n",
    "    - Преимуществом такого подхода будет в том, что нейронная сеть сама будет выявлять дополнительные признаки и зависимости, что позволит уменьшить сложность подготовительной части проекта, связанной с анализом признаков.\n",
    "    \n",
    "    "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  },
  "widgets": {
   "state": {
    "2e934aa725a34d238e9e99f51bfb9c27": {
     "views": [
      {
       "cell_index": 57
      }
     ]
    },
    "ceba48674562412188463e3a48356ee8": {
     "views": [
      {
       "cell_index": 57
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
